function [t,Tg,Tgb,ATg,N,Dcum,Nd,Ns,Ndead,r,control]=kinetic_model_Markov(param,dose0,halflife0,number_of_injections)    
%This code implements a stocastic thyroid response model with different types of cells: undamaged, sublethally damaged, lethally damaged and dead, 
%as well as two tumor markers, namely the serum concentration of Tg and TgAb. A pretreatment phase simulates the growth of tissue. This allows to
%obtain a thyroid tissue with a mass equal to some predetermined value (for example the mass of a experimentally measured thyroid/ lesion)
%and the corresponding concentration of serum Tg and TgAb. The treatment phase simulates a I-131 irradiation producing sublethal and lethal 
%damage. The treatment can involve several I-131 administrations, delivered at periodic time intervals. 
%Transfer probabilities for each compartment/process are defined for the variables corresponding to cells of different types, and used in a stochastic
%Markov model. Individual stochastic calculation is performed when the cells of each compartment are below a threshold value fixed to 10000. The 
%evolution of cell compartments when N>threshold is scaled up using the
%results from the 10000 simulated stochastic processes.

%Input: 
%1-param: vector with the model parameter values (see lines 28-48 for a description).
%2-dose0: averaged population dose achieved with a standard treatment administration.
%3-halflife0: averaged population effective halflife of the radionuclide in the tumor.
%4-number_of_injections: number of treatment administrations.

%Output:  
%1-t:        time vector
%2-Tg:       vector with the temporal evolution of free Tg serum concentration in ng/ml
%3-Tgb:      vector with the temporal evolution of Tg serum  bound to TgAb concentration in ng/ml
%4-ATg:      vector with the temporal evolution of the TgAb serum concentration in ng/ml
%5-N:        vector with the temporal evolution of the number of undamaged cells
%6-Dcum:     vector with the temporal evolution of the absorbed dose 
%7-Nd:       vector with the temporal evolution of the number of doomed cells
%8-Ns:       vector with the temporal evolution of the number of sublethally damaged cells
%9-Ndead:    vector with the temporal evolution of the number of dead cells
%10-r:       vector with the temporal evolution of the dose rate
%11-control: time before treatemtn


%Description of input parameters vector
%Proliferation parameters:
%Proliferation:
lambda=log(2)/param(1);               %Normal cells prolif rate (day^-1)
lambda_d=log(2)/param(2);             %Doomed cells prolif rate (day^-1)
lambda_s=log(2)/param(3);             %Sublethally damaged prolif rate (day^-1)
lambda_t=param(4);                    %Thyroglobulin production rate in ng/(ml·day)

%Radiosensitivity, repair & resorption rates
a=param(5);                           %Letal damage (damaged cells/decay)
b=param(6);                           %Subletal damage (damaged cells/decay)
mu=log(2)/(param(7)/24);              %Subletal damage repair rate (day^-1),(param(7) is in hours, 2h halflife)
gamma=log(2)/param(8);%=log(2)/(10);  %Resorption rate for doomed cells (Ln(2)/T_1/2, day^-1)
ke=log(2)/param(9);                   %Rate of thyroglobulin elimination (day^-1)
N0=param(10);                         %Number of undamaged cells just before treatment
cte_Tg=param(11);                     %Thyroglobuline production modulation constant (ng/ml)
Nmax=param(12);                       %Threshold for saturation of cell proliferation
gamma2=log(2)/param(13);              %Rate of death for doomed cells (Ln(2)/T_1/2, day^-1)
k=log(2)/param(14);                   %TgAb natural elimination rate (day^-1)
ka=log(2)/param(15);                  %Free Tg binding rate to Tg antibodies in ml/(IU·day), or ml/(ng day) for Tg antibodies
z=log(2)/param(16);                   %TgAb production rate in IU/(ng·day)
keb=log(2)/param(17);                 %Bound Tg enhanced metabolic clearance rate (day^-1)
zb=param(18);                         %TgAb production rate, lymphocyte related in IU/(ml·day)
zt=param(19);                         %TgAb production rate in ng/(ml·day)
halflife=param(20);                   %patient specific effective halflife of the radionuclide in the lesion (day)
uptake=param(21);                     %patient specific radionuclide uptake, this is the parameter scaled to change the administered activity 
                                      %in the patient individualization strategy. 

resting_time=180;  %Time between injections (days).


decay_rate0=log(2)/halflife0;        %averaged population decay rate of the radionuclide in the lesion
decay_rate=log(2)/halflife;          %patient specific decay rate of the radionuclide in the lesion

%Initial conditions: the simulation begins with a small number of undamaged proliferating cells and all other variables =0.

N(1)=100000;     %Initial undamaged cells
Nd(1)=0;         %Initial doomed cells
Ndead(1)=0;      %Initial dead cells
Ns(1)=0;         %Initial subletally damaged cells
Tg(1)=0;         %Initial Tg (free compartment)
Tgb(1)=0;        %Initial Tg (bound compartment)
ATg(1)=0;        %Initial TgAb
r1(1)=0;         %Initial dose rate in the interstitial compartment
r2(1)=0;         %Initial dose rate in the interstitial compartment
r(1)=r1+r2;      %Total dose rate
Dcum(1)=0;       %Initial absorbed dose
t(1)=0;          %Initial time

dt=5/(60*24);   %Time step, 20 minutes in days
i=1;

%Simulation of tissue evolution until the thyroid mass equal to the mass achieved when the treatment begins.
while N(end)<N0
     %Cell variables are updated using the MarKow model
     [N(i+1),Ns(i+1),Nd(i+1),Ndead(i+1),p_d_sub(i+1),p_d_let(i+1),p_prol1(i+1),p_prol2(i+1),p_rep(i+1),p_death(i+1),p_lysis(i+1)]=surv_Markov(N(i),Ns(i),Nd(i),Ndead(i),a,b,r(i),mu,lambda*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax),lambda_d*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax),lambda_s*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax),gamma,gamma2,dt);
    
     %Tiroglobulin (free) compartment
     dTg = -ke*Tg(i)...
         +lambda_t*(N(i)+Ns(i)+Nd(i))/(cte_Tg+(Tg(i)+Tgb(i)))...
         -ka*ATg(i)*Tg(i);

     %Tiroglobulin (bound to TgAb) compartment
     dTgb=ka*ATg(i)*Tg(i)...
         -(ke+keb)*Tgb(i);
     
    %AntiTiroglobulin autoantibodies compartment
     dATg=z*(Tg(i)+Tgb(i))...
         +zb*Nd(i)...
         -ka*ATg(i)*Tg(i)...
         -k*ATg(i);
     
    %dose rate (Equal to zero. It is defined in this growing phase to have all variable vectors with same size)
     dr  = -decay_rate*r(i);

     %Variables update
     Tg(i+1,1)=Tg(i)+dTg*dt;
     Tgb(i+1,1)=Tgb(i)+dTgb*dt;
     ATg(i+1,1)=ATg(i)+dATg*dt;
     r(i+1,1)=r(i)+dr*dt;
     Dcum(i+1)=Dcum(i)+r(i)*dt;     
     t(i+1)=t(i)+dt;
     i=i+1;
end
%Construction of the a vector with the times of treatment administration.
t_before_trat=t(end)+dt;
indice=length(t);
t_end=t_before_trat+(number_of_injections-1)*resting_time+6.67*30; 
treatment_time=[0 [t_before_trat:resting_time:t_before_trat+resting_time*(number_of_injections-1)] t_end]; 
%Variable used to control the number of treatment administrations already delivered
injection=1;
t_tot=(0:dt:t_end+dt);

%Simulation of the treatment
for j=2:(length(treatment_time)-1)

    t_span=(treatment_time(j):dt:treatment_time(j+1)-dt);
    t=[t t_span];
    steps=length(t_span);
    if injection<=number_of_injections
        %Initialization of the dose rate value at the beginning of the adminstration
        r0=dose0*decay_rate0*uptake;
        r(1,end)=r0;
    end
    
    for i=length(N):length(N)+(steps)-1
        [N(i+1),Ns(i+1),Nd(i+1),Ndead(i+1),p_d_sub(i+1),p_d_let(i+1),p_prol1(i+1),p_prol2(i+1),p_rep(i+1),p_death(i+1),p_lysis(i+1)]=surv_Markov(N(i),Ns(i),Nd(i),Ndead(i),a,b,r(i),mu,lambda*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax),lambda_d*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax),lambda_s*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax),gamma,gamma2,dt);
        %Thyroglobulin (free) compartment 
        dTg = lambda_t*(N(i)+Ns(i)+Nd(i))/(cte_Tg+(Tg(i)+Tgb(i)))...
          +zt*Nd(i)...
         -ke*Tg(i)...
         -ka*ATg(i)*Tg(i);
     
         %Thiroglobulin (bound to TgAb) compartment      
         dTgb=ka*ATg(i)*Tg(i)...
         -(ke+keb)*Tgb(i);
         %AntiThyroglobulin autoantibodies compartment
         dATg=z*(Tg(i)+Tgb(i))...
         +zb*Nd(i)...
         -ka*ATg(i)*Tg(i)...
         -k*ATg(i);     
        %dose rate    
         dr  = -decay_rate*r(i);
        
        %Variables update 
        Tg(i+1,1)=Tg(i)+dTg*dt;
        Tgb(i+1,1)=Tgb(i)+dTgb*dt;
        ATg(i+1,1)=ATg(i)+dATg*dt;
        r(i+1,1)=r(i)+dr*dt;
        Dcum(i+1)=Dcum(i)+r(i)*dt; 
    end
 %Update the administration control variable 
 injection=injection+1;
end
%Calculate if the tumor is controlled (control= 1 if the treatment killed all the undamaged cells, control=0 otherwise). 
control=N(end)./N(end);
control(find(isnan(control)))=0;
control = double(~control);
%The output corresponds to the evolution of the variables from the beginning of the treatment
t=t(indice+1:end)-t_before_trat;%+dt;
Dcum=Dcum(indice+1:end);
Tg=Tg(indice+1:end)';
Tgb=Tgb(indice+1:end)';
ATg=ATg(indice+1:end)';
N=N(indice+1:end);
Nd=Nd(indice+1:end);
Ns=Ns(indice+1:end);
Ndead=Ndead(indice+1:end);
r=r(indice+1:end)';


%Sample the vectors (1 point per day)

steps_per_day=24*60/5; %(minutes in a day)/(dt)
t=t(1:steps_per_day:end);
Tg=Tg(1:steps_per_day:end);
Tgb=Tgb(1:steps_per_day:end);
ATg=ATg(1:steps_per_day:end);
N=N(1:steps_per_day:end);
Dcum=Dcum(1:steps_per_day:end);
Nd=Nd(1:steps_per_day:end);
Ns=Ns(1:steps_per_day:end);
Ndead=Ndead(1:steps_per_day:end);
r=r(1:steps_per_day:end);

end

function [N,Ns,Nd,Ndead,p_d_sub,p_d_let,p_prol1,p_prol2,p_rep,p_death,p_lysis]=surv_Markov(N,Ns,Nd,Ndead,a,b,r0,mu,lambda,lambda2,lambda3,gamma,gamma2,dt)
n_sim=10000;
%Calculate the number of cells in each compartment, up to a maximum value of n_sim=10000
s_loop1=min(n_sim,N);
s_loop2=min(n_sim,Ns);
s_loop3=min(n_sim,Nd);
s_loop4=min(n_sim,Ndead);
%Calculate, for each compartment, the scale factor to be applied if the number of cells is larger than n_sim
factorN=max(1,N/n_sim);
factorNs=max(1,Ns/n_sim);
factorNd=max(1,Nd/n_sim);
factorNdead=max(1,Ndead/n_sim);

%Calculate the transfer probability associated to each process
p_d_sub=a*r0*dt;
p_d_let=b*r0*dt;
p_prol1=lambda*dt;
p_prol2=lambda2*dt;
p_prol2=lambda3*dt;
p_rep=mu*dt;
p_death=gamma*dt;
p_lysis=gamma2*dt;

%Calculation of the transfer of cells betwwen compartments.
%N
p_N1=p_d_sub;
p_N2=p_N1+p_d_let;
p_N3=p_N2+p_prol1;
r=rand(s_loop1,1);
%Sublet damage
%The undamaged cells for which the random number is <P_N1 go to the sublethally damaged compartment
dN=nnz(r<= p_N1); 
N=N-dN*factorN;
Ns=Ns+dN*factorN;
%Lethal damage
%The undamaged cells for which the random number is in the interval [P_N1 P_N2] go to the doomed compartment
dN = nnz(r>p_N1 & r<= p_N2);
N=N-dN*factorN;
Nd=Nd+dN*factorN;
%Proliferation
%The undamaged cells for which the random number is in the interval [P_N2 P_N3] proliferate leading to a new cell
dN = nnz(r>p_N2 & r<= p_N3);
N=N+dN*factorN;


%Ns
p_Ns1=p_prol1;
p_Ns2=p_Ns1+p_d_sub+p_d_let;
p_Ns3=p_Ns2+p_rep;
r=rand(s_loop2,1);
%Ns prolif
dNs=nnz(r<= p_Ns1); 
Ns=Ns+dNs*factorNs;
%Lethal damage to Ns
dNs = nnz(r>p_Ns1 & r<= p_Ns2);
Ns=Ns-dNs *factorNs;
Nd=Nd+dNs *factorNs;
 %repair
dNs = nnz(r>p_Ns2 & r<= p_Ns3);
Ns=Ns-dNs*factorNs;
N=N+dNs*factorNs;


%Nd
p_Nd1=p_prol2;
p_Nd2=p_Nd1+p_death;
r=rand(s_loop3,1);
%Ns prolif
dNd=nnz(r<= p_Nd1); 
Nd=Nd+dNd*factorNd;
%Death of doomed cells
dNd = nnz(r>p_Nd1 & r<= p_Nd2);
Nd=Nd-dNd*factorNd;
Ndead=Ndead+dNd*factorNd;


%Ndead
p_Ndead1=p_lysis;
r=rand(s_loop4,1);
%Ns prolif
dNdead=nnz(r<= p_Ndead1); 
Ndead=Ndead-dNdead*factorNdead;


%Round final numbers to avoid 
N=round(N);
Ns=round(Ns);
Nd=round(Nd);
Ndead=round(Ndead);


end